MODULE fefmm_driver

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_get_FF_J_matrix,      &
             fefmm_get_NF_J_matrix

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_get_FF_J_matrix(nBas,dens,J_matrix)

      USE fmm_scheme_builder,       ONLY: fmm_get_scheme
      USE fmm_driver,               ONLY: fmm_get_multipole_potential
      USE fmm_boxed_multipole_ints, ONLY: fmm_pack_boxed_mpoles,       &
                                          fefmm_contract_multipole_potential

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)    :: nBas
      REAL(REALK),   INTENT(IN)    :: dens(nBas,nBas)
      REAL(REALK),   INTENT(INOUT) :: J_matrix(nBas,nBas)

      TYPE(scheme_paras), POINTER :: scheme
      REAL(REALK), ALLOCATABLE :: Vpot(:,:)
      REAL(REALK)   :: dummy_dens(1,1) = one
      INTEGER(INTK) :: lmax, lmdim, nbox

      CALL fmm_get_scheme(scheme)
      lmax = scheme%raw_LMAX

      ! Pack moments in same box to reduce memory in FMM readin
      CALL fmm_pack_boxed_mpoles(lmax,nBas,dens,nbox)

      ! Get multipole potential for each box
      lmdim = (1 + lmax)**2
      ALLOCATE (Vpot(lmdim,nbox))
      CALL fmm_get_multipole_potential(FE_FMM,dummy_dens,Vpot)

      ! Contract potential with AO moments box by box
      CALL fefmm_contract_multipole_potential(nBas,nbox,lmax,Vpot,J_matrix)

      DEALLOCATE (Vpot)

   END SUBROUTINE fefmm_get_FF_J_matrix

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_get_NF_J_matrix(basis,dens,J_matrix)

      USE fmm_scheme_builder, ONLY: fmm_get_scheme
      USE fmm_shell_pairs,    ONLY: fmm_get_shell_pairs
      USE fefmm_box_utils,    ONLY: fefmm_get_boxes,                &
                                    fefmm_free_boxes
      USE fefmm_box_class,    ONLY: fefmm_init_box_class,           &
                                    fefmm_free_box_class,           &
                                    fefmm_update_J_matrix

      IMPLICIT NONE
      TYPE(fmm_basis), INTENT(IN)    :: basis
      REAL(REALK),     INTENT(IN)    :: dens(basis%nbas,basis%nbas)
      REAL(REALK),     INTENT(INOUT) :: J_matrix(basis%nbas,basis%nbas)

      CHARACTER(LEN=13),  PARAMETER :: Name = 'FE_AO_overlap'
      TYPE(fmm_sh_pairs), POINTER :: sh_pairs(:)
      TYPE(scheme_paras), POINTER :: scheme
      REAL(REALK),        POINTER :: boxes(:,:)
      CHARACTER(LEN=255) :: FBuf

      REAL(REALK), SAVE :: nel = zero
      REAL(REALK)   :: tmp
      INTEGER(INTK) :: i

      CALL fmm_get_scheme(scheme)
      SELECT CASE (scheme%LIPN)
         CASE (1)
            IF (MOD(scheme%FEdim-1,2)/=0) CALL fmm_quit('FEdim must be 2k+1')
         CASE (2)
            IF (MOD(scheme%FEdim-1,2)/=0) CALL fmm_quit('FEdim must be 2k+1')
         CASE (3)
            IF (MOD(scheme%FEdim-1,3)/=0) CALL fmm_quit('FEdim must be 3k+1')
         CASE (4)
            IF (MOD(scheme%FEdim-1,4)/=0) CALL fmm_quit('FEdim must be 4k+1')
         CASE (6)
            IF (MOD(scheme%FEdim-1,6)/=0) CALL fmm_quit('FEdim must be 6k+1')
         CASE DEFAULT
            CALL fmm_quit('LIP too high order!')
      END SELECT

      CALL fmm_get_shell_pairs(basis,sh_pairs)
      CALL fefmm_get_boxes(scheme%grain,sh_pairs,boxes)
      CALL fefmm_init_box_class(scheme,basis,sh_pairs,dens)

      FBuf = TRIM(Name)//".fefmm"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='OLD', ACTION='READ',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)

print *, "Computing NF J_matrix..."
      DO i = 1, SIZE(boxes,2)
!print *, "box: ", i, "of", SIZE(boxes,2)

!         CALL fefmm_build_potential(boxes(:,i),tmp)
!         IF (ABS(tmp) < 1e-15) CYCLE   !! broken????
!         CALL fefmm_contract_FE_potential(scheme,basis,sh_pairs,    &
!                                          boxes(:,i),J_matrix)
         CALL fefmm_update_J_matrix(boxes(:,i),tmp,J_matrix,basis%nbas)
         nel = nel + tmp

      END DO
print *, "Total # electrons  = ", nel

      NULLIFY (boxes)
      CLOSE(UNIT=IO2, STATUS='KEEP')
      CALL fefmm_free_box_class
      CALL fefmm_free_boxes

   END SUBROUTINE fefmm_get_NF_J_matrix

!-------------------------------------------------------------------------------

END MODULE fefmm_driver
